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We present a rigorous and efficient approach to the calculation of classical lattice-dynamical 
quantities from simulations that do not require an explicit solution of the time evolution. We focus 
on the temperature-dependent vibrational spectrum. We start from the moment expansion of the 
relevant time-correlation function for a many-body system, and show that it can be conveniently 
rewritten by using a basis in which the low-order moments are diagonal. This allows us to compute 
the main spectral features (e.g., position and width of the phonon peaks) from thermal averages 
available from any statistical simulation. We successfully apply our method to a model system that 
presents a structural transition and strongly temperature-dependent phonons. Our theory clarifies 
the status of previous heuristic schemes to estimate phonon frequencies. 



PACS numbers: 63.10.-|-a, 63.70.-fh, 63.20.dk 

Is it possible to compute any equilibrium property of 
a material from appropriate thermal averages? From 
an atomistic simulation perspective, an affirmative an- 
swer to this question implies that, if we can characterize 
the configuration space accessible at a temperature T 
- as can be done e.g. with Monte Carlo (MC) meth- 
ods -, we can also calculate any quantity of interest. 
In particular, time-correlation functions would become 
available, which would allow us to compute T-dependent 
phonon spectra or lattice thermal transport coefficients, 
to name two non-trivial properties, without explicitly 
solving the equations of motion. The advantages of such 
a non- dynamical approach would be numerous: as in MC 
schemes, fixing T would become trivial, and we would 
avoid the problems associated with the use of thermostats 
in molecular dynamics (MD) simulations; there would be 
no need for very long MD runs to access low-frequency 
phenomena, etc. 

In the particular case of the vibrational spectrum, sev- 
eral heuristic schemes have been proposed to realize such 
a goal. For example, Hellman et al. |1| compute phonons 
at finite T from an effective dynamical matrix that they 
postulate and obtain in a computationally efficient way; 
a similar effective potential is used in Refs. i and i to 
capture anharmonic effects in the vibrational spectrum. 
While such methods are valuable, the status of the ap- 
proximations involved is not clear, and we have to resort 
to earlier works to find a more rigorous treatment. In- 
deed, as recognized in a variety of fields a time- 
correlation function can be obtained from the knowledge 
of its moments, which are static quantities that can be 
computed from efficient statistical simulations. In the 
context of lattice-dynamical studies, the moment-based 
approach introduced by Mori 0] has been used to in- 
vestigate a number of simple systems at both the quan- 
tum and classical levels |7Hll|: however, it has failed to 



gain popularity in the field of classical MD. As far as 
we can see, this is probably because (1) we lack a gen- 
eral scheme to tackle arbitrarily complex materials with 
many atoms in the cell, (2) there has not been enough 



work to identify the simplest approximations that may 
render useful results (i.e., reliable phonon frequencies and 
possibly peak widths), and (3) these ideas have remained 
much confined within the quantum-dynamics community. 
Here we present our own derivation of a moment-based 
formalism for the classical case, remedying the mentioned 
deficiencies. The resulting theory allows us to clarify the 
status of the heuristic methods in the literature [ij, 0] , 
and in our opinion should become a standard tool in the 
field of classical simulations. 

General formalism.- We define the classical correlation 
function of quantities Ai and Bj as 



C^^{t) = {MO)B,{t)) 



C^^(a.)e-*d^ (1) 



where (...) indicates thermal averaging and C^^{uj) is 
the corresponding spectral function. Below we will iden- 
tify Ai and Bj with atomic positions or velocities, i and 
J being composite indices that label an atom and a di- 
rection in space. Note that C:^^{t) is real, which implies 



Let us define C:^{t) 



For simplicity, we 



will work with the real part C{j{ijo) = 9fi[Cj^(aj)], as this 
will contain the information about the vibrational spec- 
trum Thus, in the time domain we have 
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which is even with respect to time inversion and can be 
Taylor expanded in the fohowing way 
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with the moments given by 



(2n!) 



A,2n 
My 



w""C,^(a;)da;. 



(3) 



(4) 



2 



The next key step is to prove that these moments can 
be written as certain correlation functions at t = 0, and 
can thus be computed as regular thermal averages. By 
successive partial integration, one can see that 
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we have in writing the moments, note that alternatively 
we can get 



X .4 



'(3) 



(10) 



It is also interesting to note that other time-correlation 
functions can be readily computed from Eq. 0. Indeed, 
it can be seen from Eqs. (|3]) and ([5]) that 



where A]^' is the n-th time derivative of Ai. These iden- 
tities allow us to rewrite Eq. Q in different ways. From 
a computational viewpoint, it is convenient to use 
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which involves time derivatives of the lowest possible 
order. Finally, the derivatives can be computed using 
Hamilton's equation of motion 



{A m - V _ dAdH_ 

^ \dxi dpi dpi dxi 



(7) 



The Hamiltonian is 



(8) 



where m^, Xi, and pi are mass, position and momentum, 
respectively, and F is a velocity-independent potential. 
It is convenient to introduce x[ = ^/mlxi, to get rid of 
the mass dependence in the kinetic energy. We will thus 
work with the moments /i^- of the correlation func- 



tion Cf!,'{t) 



— yfmimjC^{t), where A'^ = ^JralA^. The 
simplest object that gives information about the vibra- 
tional spectrum is the position-position correlation func- 
tion. By taking A!^ = x\ = x\ — (x'J, we obtain 



= y^milTlj {{XiXj) - {Xi){Xj)) , 
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for the four lowest moments of Eq. ^ . Here we reverted 
the mass scaling at the last step of the derivation, so as 
to express the moments in terms of the regular atomic 
positions. We also used (p-p^) = k^TSij = (i~^5ij and 
{fm})9{{<})) = (/(W}))(5({a;a)): where/and^are 
arbitrary functions. Let us stress that this procedure ren- 
ders /^fj'^ proportional to the identity matrix, a fact that 
will be advantageous later. As an example of the freedom 



,2n 



A,2n+2 



(11) 



Thus, for example, the lowest non-diagonal moment of 
the velocity- velocity correlation function is /i"^'^ = Mij '^- 
In order to gain physical insight, consider the case of 
a single harmonic oscillator. There, we know the exact 
form of (w) = /i^ ''^[S{uj~ujo) + 6{oj + ujo)]/2, and from 
Eqs. dH) and (O we obtain Wq" = then, for 

n = 1 we have mujQ — [I3{{x'^) ~ (x)^)]^-^. We thus find 
that the fluctuations of the position give us the frequency 
of the (phonon) peak in the (w) spectrum. Similarly, 
for iuj) we get mwQ = {d^V/dx^), i.e., the frequency is 
given by the thermal-averaged dynamical matrix in this 
case. These intuitive relations, which are exact in the 
harmonic limit and are generalized below to the many- 
body case, are implicitly underlying the heuristic meth- 
ods to estimate phonon frequencies mentioned above [H- 

Practical scheme - In general we will have many inter- 
acting atoms and an anharmonic potential V . To simplify 
the problem, let us make a unitary coordinate transfor- 
mation A'^ — T*^A^ that will be analogous to the 
usual change into a normal-mode basis. We write the 
transformed correlation function and moments as 
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We choose Tai to diagonalize the lowest non-diagonal mo- 



ment fif, , so that 
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Note that Tai gives the polarization vectors of our normal 
modes. We then approximate [13] 



TLC^bi^b 



ah 



TLC^l{t)T,a, (15) 



where the second equality is exact to low order in the 
moment expansion. Hence, to investigate the spectrum 
given by {uj), we will work with the collection of an- 
harmonic oscillators C^a (w) [ijj 
Now, we have jlfj 



and fil^ 



/i''/^. Interest- 
ingly, thanks to the mass-scaling transformation, fxi^' 
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and are proportional to the identity matrix and will 
remain diagonal in our normal-mode basis. Thus, since 
the two lowest-order moments are diagonal for both Cf^ 

and C^j , it seems reasonable to propose the following 
effective harmonic approximation 
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A', 2 
f^aa 



x',0 



I A'.O 
/^J.aa ■ 

Eq. (fT4)) involves the diagonalization of 



y/mimj{{xiXj) — {xi){xj)); for A[ — v[ we di- 

agonalize Mij '^ = which can be expressed as the 

thermal-averaged dynamical matrix of Eq. (jlOl) . Hence, 
the eigenvalues of these matrices provide us with a rig- 
orously justified approximation to the position of the 
(phonon) peaks in the C^^ {uj) spectra; as described 
above, they render the exact phonon frequencies in the 
harmonic limit. 

In order to capture more complex line shapes, we need 
a strategy to treat the full C^(a;). One possibility is to 
assume an analytic form for it, with free parameters that 
can typically be written as a function of the low-order 
moments. We worked with some physically motivated 
choices, i.e., Gaussians, Lorentzians, and combinations of 
delta functions. Such an approach allowed us to compute 
the peak frequencies for a variety of model systems, even 
in the presence of significant non-linear effects (i.e., over- 
tones). However, the scheme failed to provide a quanti- 
tative estimate of the peak widths; further, for models 
with very strongly T-dependent frequencies, we some- 
times obtained non-physical solutions for the parameters 
of the trial spectral functions. 

Fortunately, we found it possible to treat C^^j {uj) 
in a robust and accurate way by resorting to Mori's 
continued- fraction representation [6| . Using the notation 
ofRef.i, we write for a generic, real autocorrelation func- 
tion C{oj) 
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where the (5„ parameters are explicit functions of the mo- 
ments. The first three are given by Q 
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A continued fraction is usually terminated by assuming 
the last ipn{z) term to be the Laplace transform of some 
model function. Many terminations have been used in 
the literature [7l-[ll|: yet, we found that, for the model 



systems we investigated, the line shape is quite insen- 
sitive to the termination scheme when including up to 
sixth moments; thus, we simply used for -02(2) the Gaus- 
sian termination described in Ref. S Obviously, in this 
scheme a separate continued-fraction expansion must be 
calculated for each C^(a;). In practice, we first com- 
pute the atomic moments /i^ and then obtain those 
corresponding to our normal modes by using Eq. ()13|) . 

Example of application.- To test our approach, we used 
it to compute vibrational spectra - with moments calcu- 
lated from MC simulations - and compared the results 
with the exact ones obtained from MD 14, [la. We 



worked with model systems, which gave us full control 
of the potential-energy surface and allowed us to try the 
method in very diverse and challenging situations. Over- 
all we found that our approach renders excellent results 
for the main features of the vibrational spectrum. Here 
we describe a particularly demanding case, namely, a sys- 
tem undergoing a structural phase transition driven by a 
soft phonon mode. 

Let us consider a cubic crystal with three degrees of 
freedom xia per cell I, where a — x, y, and z. For sim- 
plicity we take m = 1 and write the Hamiltonian as 



E/2 2, 2 2, 2 2\, 
y^lx^ly + ^ly^lz + ^Iz^lx) + 



(19) 



a' 



Calf; - X;' P + C4 ^ \xia - Xi'J 



where the primed sum is restricted to nearest-neighboring 
cells. In essence, this is the well-known discrete (j)^ model 
[l6| . with the following extensions: we included an on- 
site anisotropic term (c2 > 0) chosen so that the ground 
state has a tetragonal symmetry, and a coupling between 
nearest-neighbors {cf) that breaks the symmetry between 
longitudinal and transversal phonon branches. Here we 
show representative results obtained for a choice of pa- 
rameters (ci = 0.25, C2 = 0.50, C3 = 1.00, a = 0.50) that 
leads to a second-order displacive [l^l phase transition. 
While the energy units are arbitrary, this model renders 
a realistic representation of a phase transition at room 
temperature. 

We simulated the model in a periodically-repeated 
20x20x20 simulation box. We carefully checked the con- 
vergence of the MC and MD simulations. For example, 
the MD results shown here were obtained by Fourier 
transforming time-correlation functions computed from 
constant-energy trajectories whose starting points were 
snapshots taken from a constant-T Langevin simulation. 
For each T investigated, ten different starting points were 
considered, the final spectral function being an average. 

In periodic systems the moments, as well as the 
time-correlation and spectral functions, become block- 
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FIG. 1. (Color online.) Dynamical properties of our model 
system. We show the results obtained from our approach 
based on IVIC simulations (blue) and the exact MD results 
(red). Panels (a) and (b): Phonon bands at two different T's. 
Panel (c): Spectrum at the F point, for T's around the critical 
temperature Tq. Panel (d): F spectrum at selected T's. All 
spectra are normalized to unity. See text for more details. 



diagonal in the Bloch representation. Hence, when dis- 



cussing our model, we can substitute C. 



A' 



C. 



A' 

qQ,-q0 



c. 



c. 



q.aa ' 



etc., where q is 



within the first Brillouin zone. 

Figures lUa) and (Hb) show the phonon bands at low 
and high T's, respectively. The MD results were obtained 
from the first moment of the peaks in the spectrum given 
by the trace of ^^^(w). The JMC results were obtained 

by diagonalizing jj"^ The agreement is excellent. 

Figure [TJc) shows the T-dependent frequencies of the 
phonons at the F point (q = 0). In the cubic phase, 
the frequencies decrease as T is reduced, and essentially 
vanish at the critical temperature Tq. Then, below Tc 
the frequencies increase as T decreases. We have three F 
phonons in both phases: they are three-fold degenerate 
in the cubic structure, but split in two groups when the 
symmetry is lowered to tetragonal. In the figure we show 
the trace of the Cf, ^^{uj) functions resulting from MD 
simulations, as well as the frequencies computed by diag- 



onalizing /ip obtained from MC. This approximation 
gives excellent results, even in the immediate vicinity of 
Tc where the system is strongly anharmonic. 

Figure (Hd) shows the line shapes from the continued- 
fraction representation of Cp aai'^)^ using up to the sixth- 
order moments, together with the MD results. We can 
appreciate that the widths of the peaks obtained with 
our method are semi-quantiatively correct. We include 
a result at an unrealistically high T = 100, where the 
peak broadening is very significant. Even in such extreme 
conditions, our approximate spectral function provides a 
fair representation of the exact one. 

Final remarks - We have shown that the main fea- 
tures of classical vibrational spectra can be accurately 
computed from knowledge of the low-order moments of 
the appropriate time-correlation functions. The moments 
can be obtained as thermal averages from MC simula- 
tions. Alternatively, one may obtain them from MD 
simulations, without the need to explicitly compute the 
time-correlation functions; this should allow for shorter 
MD runs (only as long as needed to compute accurate 
thermal averages) and simplify the use of thermostats 
(as their interfering with the dynamics would be unim- 
portant). 

The present method renders accurate results for the 
vibrational frequencies. Remarkably, very good results 
are already obtained within the crudest approximation, 
which amounts to considering the effective harmonic sys- 
tem defined by Eq. (fTC)) . Indeed, our method allows us to 
obtain accurate vibrational frequencies from appropriate 
thermal averages of atomic positions or forces, which are 
readily available in any atomistic simulation. Further, 
our effective-harmonic treatment provides a rigorous jus- 
tification for some of the assumptions underlying previ- 
ous schemes in the literature We have also shown 
that it is possible to reproduce the line shape of the spec- 
tral functions in a semi-quantitative way, provided higher 
moments are available. 

We hope the methods here discussed will become stan- 
dard tools in classical simulations, where they can be 
used to a great advantage. Further, we hope this work 
will serve as an example of the usefulness of a formalism 
that, unfortunately, seems to remain unknown to many in 
the MD community. Consider, for example, the identity 
{d'^V/dxidxj) = (3{{dV/dxi){dV/dxj)), which we have 
proven here and constitutes a rigorous way to compute 
the thermal-averaged force-constant matrix from appro- 
priate products of forces. Such a relation should be very 
useful in the context of first-principles simulations, where 
the calculation of second derivatives of the energy may 
be computationally prohibitive. For instance, it would 
simplify, and provide a solid foundation for, the practical 
implementation of the scheme of Hellman et al. [if, to 
name one obvious possibility. 
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